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1.  Introduction 


Synthetic  2-phase  microstructures  are  often  used  as  a  surrogate  for  real  microstruc¬ 
tures  based  on  experimental  microstructure  statistics  or  as  a  way  to  test  how  dif¬ 
ferent  micro  structure  features  affect  properties  of  materials.  These  can  be  gener¬ 
ated  through  processes  such  as  the  random  sequential  algorithm  described  within 
or  other  techniques,  some  of  which  are  outlined  in  the  book  on  heterogeneous  mi¬ 
crostructures  by  Torquato.1  Random  sequential  adsorption  (RSA)  has  been  applied 
to  understand  microstructure  and  phenomena  in  a  number  of  different  systems  and 
applications,  including  far-from-equilibrium  processes,2  car  parking  and  protein  ad¬ 
sorption,3  colloids,4  dimer  adsorption,5  fiber  composites,6  and  even  disks,7  rectan¬ 
gles8-9  and  n-star  objects.10 

In  this  technical  note,  a  random  sequential  adsorption1112  algorithm  implemented  in 
MATLAB13  is  used  to  generate  three-dimensional  (3-D)  synthetic  microstructures 
composed  of  packing  of  ellipses  within  a  voxelized  microstructure.  The  MATLAB 
scripts  are  attached  as  appendices  for  future  development  and  application.  This 
work  is  an  extension  of  a  previous  two-dimensional  (2-D)  synthetic  microstructure 
builder  that  incorporates  2-D  ellipses,  aspect  ratio,  area  fraction,  size  (polydisper- 
sivity),  orientation,  and  size/orientation  distributions.14 

2.  Algorithm/Code  Description 

The  algorithm  used  to  generate  the  synthetic  microstructures  is  an  RSA  algorithm, 
which  is  typically  used  to  represent  particle  deposition  on  surfaces.  For  instance, 
there  are  multiple  stages  in  particle  deposition:  1)  initially,  the  surface  is  clean  and 
free  of  particles,  2)  then  particles  begin  to  spontaneously  attach  to  the  surface  at 
small  concentrations,  and  3)  the  deposition  progressively  slows  down  due  to  de¬ 
posited  particles  blocking  subsequent  deposition  of  other  particles.  This  process  is 
typically  studied  by  the  RSA  model  (for  more  information,  see  Torquato1),  which 
occurs  via  the  following  steps: 

•  A  spherical/ellipsoidal  particle  is  randomly  placed  on  the  surface.  Once  placed, 
this  particle’s  position  is  fixed  and  it  cannot  move. 

•  A  trial  spherical/ellipsoidal  particle  is  randomly  placed.  If  this  particle  over¬ 
laps  previously  placed  particle(s),  then  this  attempt  is  rejected.  Otherwise,  the 
placement  is  accepted. 

Approved  for  public  release;  distribution  is  unlimited. 


1 


•  The  model  generally  proceeds  by  adding  until  a  predetermined  area/volume 
fraction  of  particles  is  reached  or  the  jamming  limit  is  reached.  The  jamming 
limit  or  saturation  of  the  surface  occurs  when  it  is  no  longer  possible  to  place 
a  particle  within  the  holes  of  the  other  deposited  particles. 

3.  Implementation  and  Usage 

This  is  implemented  through  the  MATLAB  scripts  in  Appendix  A,  Appendix  B, 
and  Appendix  C  by  using  3-D  matrices,  where  the  background  is  0  and  the  particle 
is  1.  For  the  3-D  ellipses,  it  is  possible  to  alter  the  3  orthogonal  axes  of  the  ellipse 
(i.e.,  spheres,  plates,  or  needles),  the  orientations  of  the  ellipses,  and  the  volume 
fraction.  In  both  cases,  it  is  relatively  straightforward  to  include  distributions  of 
ellipse  shapes,  ellipse  orientations,  or  ellipse  sizes.  Appendix  A  is  the  main  code, 
Appendix  B  is  the  optimized  ellipse  builder  function,  and  Appendix  C  includes  a 
few  functions  called  within  the  script. 

The  attached  scripts  have  been  tested  on  MATLAB  R2014  and  R2015  on  a  Win¬ 
dows  operating  system.  The  code  can  be  executed  by  following  these  steps: 

•  Download  the  various  scripts  into  the  same  directory: 

-  generate_synthetic_structure_ellipse_3D_anisotropic.m  (download) 

-  image_ellipse_3D_fast.m  (download) 

-  image_view.m  (download) 

-  deleteEmptyExcelSheets.m  (download) 

•  Open  the  script  “generate_synthetic_structure_ellipse_3D_anisotropic.m”  in 
MATLAB 

•  Type  “generate_synthetic_structure_ellipse_3D_anisotropic”  at  the  command 
prompt  to  run. 

Different  3-D  micro  structures  can  be  obtained  by  adjusting  the  initialization  param¬ 
eters  volume  fraction  (Vfjnax),  ellipse  axis  lengths  (a,b,c),  and  image  size  (e.g., 
the  “I=false(256,256,256);”  statement).  The  approximate  number  of  particles  are 
calculated  to  give  an  approximate  indication  of  the  time  that  may  be  required.  For 
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higher- volume  fractions,  the  code  will  have  a  progressively  harder  time  placing  par¬ 
ticles  as  less  matrix  space  is  available.  The  code  supplied  can  be  easily  extended  to 
incorporate  size  distributions,  aspect  ratio  distributions,  or  orientation  distributions, 
as  has  been  done  in  previous  work.14 

4.  Examples 

This  code  was  developed  to  examine  the  structure-property  relationship  based  on 
the  influence  of  particle  aspect  ratio,  area  fraction,  packing/clustering,  and  parti¬ 
cle  orientation.  For  instance,  prior  work  has  examined  metrics  for  characterizing 
length  scales  in  different  composite  microstructures.1516  The  voxelized  microstruc¬ 
tures  can  be  output  as  images  or  meshes  for  subsequent  calculations.  Also,  the  pa¬ 
rameters  used  to  generate  the  synthetic  microstructures  can  be  output  for  building 
within  another  software  package  (i.e.,  box  dimensions  in  x,y,z  coordinates,  center 
of  ellipse  in  x,y,z  coordinates,  orthogonal  major/minor  axis  lengths  of  the  ellipse 
[aspect  ratio],  and  the  orientation  of  the  ellipse  in  Euler  angles  F°r  in¬ 

stance,  Fig.  1  shows  an  example  of  a  subset  of  2-D  slices  from  a  3-D  synthetic  mi¬ 
crostructure  with  over  13,000  particles  placed  at  a  volume  fraction  of  30%,  where 
the  aspect  ratio  is  increased  from  spherical  (1 : 1 : 1)  to  more  needle-like  (4: 1:1).  Since 
all  ellipses  in  these  microstructures  are  identical,  Fig.  lc  shows  how  difficult  it  is  to 
deduce  the  3-D  microstructure  from  2-D  microstructure  slices  (i.e.,  this  microstruc¬ 
ture  appears  to  have  particles  with  a  distribution  of  aspect  ratios).  Figure  2  illustrates 
the  full  2-D  slice  for  this  microstructure  to  illustrate  this  more  clearly.  Notice  that 
the  periodic  boundaries  are  somewhat  visible  in  this  image  for  some  of  the  particles 
where  the  major  axis  of  the  ellipse  is  perpendicular  to  the  viewing  direction  (i.e., 
they  appear  more  elliptical  in  the  cross-section  image).  Figure  3  illustrates  the  full 
2-D  slice  for  a  microstructure  with  randomly  oriented  ellipses  that  have  an  8:8:1 
aspect  ratio  (i.e.,  plates).  Again,  the  2-D  image  appears  to  contain  a  mix  of  particle 
sizes  and  shapes  when  in  reality,  these  particles  are  all  the  same  size. 
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Fig.  1  2-D  images  of  3-D  synthetic  microstructures  with  uniform  size  circles/ellipses  that  are 
randomly  oriented  with  aspect  ratios  of  a)  1:1:1,  b)  2:1:1,  and  c)  4:1:1.  The  volume  fraction  in 
this  microstructure  is  30%  for  all  microstructures.  For  the  volume  computed,  there  were  over 
13,000  particles  placed  for  each  of  the  cases. 


Approved  for  public  release;  distribution  is  unlimited. 


4 


tfjiv 


Fig.  3  2-D  images  of  3-D  periodic  synthetic  microstructure  with  uniform  ellipses  of  8:8:1  as¬ 
pect  ratio  (i.e.,  plates)  that  are  randomly  oriented.  The  volume  fraction  in  this  microstructure 
is  19%.  For  the  volume  computed,  there  were  1,326  particles  placed. 


Approved  for  public  release;  distribution  is  unlimited. 


6 


5.  Conclusion 


Synthetic  2-phase  microstructures  are  often  used  as  a  surrogate  for  real  microstruc¬ 
tures  based  on  experimental  microstructure  statistics  or  as  a  way  to  test  how  dif¬ 
ferent  micro  structure  features  affect  properties  of  materials.  In  this  technical  note, 
an  RSA  algorithm  was  implemented  in  MATLAB  to  generate  3-D  synthetic  mi¬ 
crostructures  composed  of  packing  of  ellipses  within  a  voxelized  microstructure. 
The  MATLAB  scripts  are  attached  as  appendices  for  future  application  or  develop¬ 
ment. 
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Appendix  A.  MATLAB  Scripts:  Main  Code 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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%%  Generate  Synthetic  3D  Structure  with  Ellipses 
%  Mark  Tschopp 
%  2016 

%  Random  Sequential  Adsorption  algorithm  that  varies 
%  ellipse  volume  fraction,  size,  aspect  ratio,  and 
%  orientation 

clear  all;  clc; 

%  Initialize  variables  used  in  generation  of  3D  microstructure 
tic 

Vf  =  0;  Vf_max  =  0.05; 
a  =  1 2 ;  b=6;  c=6; 
nparticles  =  0; 

I  =  false (256,  256,  256)  ; 
image_size  =  size  (I); 

%  Calculate  approximate  Vf  of  particle  and  approximate  number  of 
%  particles  to  generate  the  specified  volume  fraction.  This 
%  will  impact  the  time  that  it  takes  to  render  the  3D  volume. 

%  Assign  random  orientation 

psil  =  2*pi*rand;  psi2  =  2*pi*rand;  phi  =  acos (rand) ; 

I_ellipse  =  image_ellipse_3D_f ast (a, b, c, psil , psi2 , phi ) ; 

Vf_particle  =  sum ( I_ellipse ( : ) )  /  numel(I); 

approx_nparticles  =  Vf_max/Vf_particle; 

disp ([' Approximate  #  of  particles:  '  num2str (approx_nparticles) ] ) ; 

%  Main  loop  for  generating  microstructure 
h  =  waitbar ( 0 ,' Ellipsoid  placement'); 
while  Vf  <  Vf_max 

psil  =  2*pi*rand; 
psi2  =  2*pi*rand; 
phi  =  acos (rand); 

I_ellipse  =  image_ellipse_3D_f ast (a, b, c, psil , psi2 , phi ) ; 

xO  =  ceil (rand*image_size ( 1 ) ) ; 
yO  =  ceil (rand*image_size (2 ) ) ; 
zO  =  ceil (rand*image_size (3) ) ; 

diam  =  size (I_ellipse, 1) ; 
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%  Merge  circle  matrix  with  synthetic  microstructure  image 
nlo  =  xO  -  floor (diam/2) ;  nhi  =  xO  +  ceil (diam/2) -1; 
ix  =  mod (nlo : nhi , image_size ( 1 ) ) ; ix ( ix==0 ) =image_size ( 1 ) ; 
nlo  =  yO  -  floor (diam/2) ;  nhi  =  yO  +  ceil (diam/2 ) -1 ; 
iy  =  mod (nlo : nhi , image_size (2 ) ) ; iy ( iy==0 ) =image_size (2 )  ; 
nlo  =  zO  -  floor (diam/2) ;  nhi  =  zO  +  ceil (diam/2 ) -1 ; 
iz  =  mod (nlo : nhi , image_size ( 3 ) ) ; iz ( iz==0 ) =image_size (3 ) ; 

Itest  =  logical (I (ix, iy, iz) ) ; 
if  sum ( Itest ( I_ellipse ) )  ==  0 
Itest ( I_ellipse )  =  1; 

I(ix,iy,iz)  =  Itest; 
accept  =  1; 

else 

accept  =  0; 
iter  =  0; 

while  accept  ==  0  &&  iter  <=  10 
iter  =  iter  +  1; 
xO  =  ceil (rand*image_size (1) ) ; 
yO  =  ceil (rand*image_size (2 ) )  ; 
zO  =  ceil (rand*image_size (3) )  ; 
diam  =  size (I_ellipse, 1) ; 


%  Merge  circle  matrix  with  synthetic 

nlo  =  xO  -  floor (diam/2 ) ;  nhi  =  xO  + 

ix  =  mod (nlo : nhi , image_size ( 1 ) )  ; 

ix  (ix==0) =image_size (1)  ; 

nlo  =  yO  -  floor (diam/2 ) ;  nhi  =  yO  + 

iy  =  mod (nlo : nhi , image_size ( 2 ) ) ; 

iy (iy==0) =image_size (2)  ; 

nlo  =  zO  -  floor (diam/2 ) ;  nhi  =  zO  + 

iz  =  mod (nlo : nhi, image_size (3) ) ; 

iz  (iz==0) =image_size (3)  ; 


micro structure 
ceil (diam/2 ) -1 ; 


ceil (diam/2 ) -1 ; 


ceil (diam/2 ) -1 ; 


Itest  =  logical ( I (ix, iy, iz) ) ; 
if  sum ( Itest ( I_ellipse ) )  ==  0 
Itest (I_ellipse)  =  1; 
I(ix,iy,iz)  =  Itest; 
accept  =  1; 

end 

end 

end 
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/  numel ( I ) ; 


if  accept; 

nparticles  =  nparticles  +  1; 

Vf_particle  =  sum ( I_ellipse ( : ) ) 

Vf  =  Vf  +  Vf_particle; 
ellipse (nparticles, 1 : 10)  =  ... 

[xO  yO  zO  a  b  c  psil  psi2  phi  Vf_particle] ; 

end 

waitbar (Vf /Vf_max) 

end 

close (h) 

disp(['3D  Digital  Slices  (sec):  '  num2str (toe)  ]  )  ; 
disp ([' Number  of  particles:  '  num2str (nparticles  )])  ; 

%%  Write  ellipse  parameters  to  Excel 

%  Store  all  ellipse  parameters  in  Excel  so  that  the  3D  structure 
%  can  be  generated  without  worrying  about  overlap. 

image_save_f lag  =  1; 
if  image_save_f lag 

string_path  =  pwd; 

Excel_f ileName  =  sprintf ( ' %s\\3D_ellipses . xlsx ' , string_path) ; 
sheet  =  sprintf (' ellipse_%dvf_%da_%db_%dc_%05d ',.. . 
round (100 *Vf ) ,a,b,c, nparticles ) ; 
warning  off  MATLAB : xlswrite :AddSheet ; 

Col Headers  =  { ' xO 1 ,  ' yO ' ,  ' z  0 ' ,  'a',  ' b ' ,  ' c ' ,  'psil',  'psil',  'phi' }; 
xlswrite (Excel_f ileName,  ColHeaders,  sheet,  'Al'); 
xlswrite (Excel_f ileName,  ellipse,  sheet,  'A2'); 
deleteEmptyExcelSheets (Excel_f ileName) ; 


end 

%%  View  slices  of  a  3D  section 

%  This  routine  views  2D  orthogonal  planes  in  the  3D  structure  as  a 
%  function  of  depth 

j=i; 

close  all; 
figure ; 
while  j 

for  i  =  1 : image_size (3 ) 

J  =  1-1  (  :  ,  :  , i) ; 
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image_view ( J)  ; 

title ([' Slice  '  num2str(i)]) 
pause  (0.05) 

end 

end 

%  Hit  control-C  to  exit 

%%  Save  individual  slices  as  jpegs  for  animated  gifs 

%  For  use  with  ImageJ  -  First,  open  ImageJ  and  select  Plugins  ->  List 
%  Opener,  then  select  the  *.txt  file  that  contains  all  the  image 
%  paths .  Once  the  images  are  loaded  in  ImageJ,  then  convert  them  to 
%  Slices  using  the  Image  ->  Stacks  ->  Image  to  Stacks  command.  Now 
%  the  stack  of  images  can  be  saved  as  an  animated  *.gif  format  for 
%  insertion  into  Powerpoint.  The  Volume  viewer  and  Volume J  renderers 
%  in  ImageJ  also  accept  stacks  of  images. 

tic 

close  all;  figure; 

sheetname  =  sprintf ( ' ellipse_%dvf_%da_%db_%dc_%05d ' , . . . 

round (100 *Vf ) , a, b, c, nparticles )  ; 
if  ~isdir (sheetname) ,  mkdir ( sheetname ) ;  end 
copyfile ( ' image_view.m' , sheetname) 
cd ( sheetname ) 

fid  =  f open ([ sheetname  ' . txt ' ] , ' w ' ) ; 
for  i  =  1 : image_size (3) 

K  =  reshape (1-1 ( : , : , i) , image_size (1) , image_size (2) ) ; 

J  =  K; 

image_view ( J) ; 
pause (0.05) 

image_f ilename  =  [sheetname,  sprintf ( '_%03d. jpg' ,  i) ] ; 
fprintf(fid,  '%s  \n ', image_f ilename)  ; 

if  image_save_f lag;  imwrite (J, image_f ilename, ' jpg ') ;  end 

end 

f close (fid) ; 

delete ( ' image_view.m' ) 

cd  .  . 

%  Save  to  *.mat  file  as  well 
save ( [sheetname  ' .mat ' ] , ' I ' ) 

disp ([' Image  viewing  and  saving  (sec):  '  num2str (toe) ] ) ; 
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%  *.mat  file  is  useful  for  subsequent  operations  in  MATLAB  or  writing 
%  files,  * . jpg  is  useful  for  visualization,  *.tif  binaries  are  a  useful 
%  if  lossless  compression  is  desired,  *.txt  contains  a  list  of  images, 

%  *.xls  contains  spreadsheet  with  all  ellipse  parameters,  if  desired 
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Appendix  B.  MATLAB  Function:  Create  Ellipse 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  y  =  image_ellipse_3D_f ast (a,  b,  c,  psil,  psi2,  phi) 

%  Mark  Tschopp 
%  2016 

%  This  subroutine  generates  a  3D  ellipse  in  a  matrix  (y)  using  the 
%  orthogonal  distances  in  the  ellipse  (a,b,c)  and  the  Euler  angles 
%  (psil , psi2, phi) ,  which  define  the  rotation  of  the  ellipse. 

%{ 

%  These  parameters  can  be  used  for  debug  purposes . 
a  =  40; 
b  =  10; 
c  =  10; 

a  =  ceil (rand*40) ; 
b  =  ceil (rand*40) ; 
c  =  ceil (rand*40) ; 
psil  =  2*pi*rand; 
psi2  =  2*pi*rand; 
phi  =  acos (rand) ; 

%} 

%  First,  identify  what  the  size  of  the  matrix  must  be  to  contain  the 
%  ellipse  and  also  what  the  central  symmetry  plane  will  be  (ic) . 

if  a  >  b;  diameter  =  2*a;  else  diameter  =  2*b;  end 
if  2*c  >  diameter;  diameter  =  2*c;  end 

if  mod (diameter, round (diameter) ) ~=  0;  diameter  =  ceil (diameter) ;  end 

radius  =  diameter/2; 

if  mod ( radius , round ( radius )) ~=  0; 

diameter  =  diameter+1;  radius  =  diameter/2; 

end 

dist ( 1 : diameter+1 , 1 : diameter+1 , 1 : diameter+1 )  =  2; 
ic  =  1  +  radius; 

%  M  is  the  rotation  matrix  for  the  ellipse,  based  on  the  Euler  angles 
%  psil,  psi2,  and  phi.  Different  rotation  matrices  and  Euler  angle 
%  conventions  can  be  used,  if  necessary. 

M  =  zeros  ( 3 ) ; 

M(l,l)  =  cos (psil)  *cos (psi2) -sin (psil) *sin (psi2) *cos (phi) ; 

M(2,l)  =  -cos (psil) *sin (psi2) -sin (psil) *cos (psi2) *cos (phi)  ; 

M(3,l)  =  sin  (psil ) *sin (phi ) ; 

M(l,2)  =  sin  (psil) *cos (psi2) +cos (psil) *sin (psi2) *cos (phi) ; 
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M(2,2)  =  -sin  (psil) *sin (psi2) +cos (psil) *cos (psi2) *cos (phi) ; 

M(3,2)  =  -cos (psil ) *sin (phi) ; 

M(l,3)  =  sin  (psi2) *sin (phi) ; 

M(2,3)  =  cos  (psi2 ) *sin (phi )  ; 

M  (3, 3)  =  cos (phi) ; 

%  This  step  creates  a  matrix  with  the  distances  to  the  ellipse 
%  centroid  (ic,  ic,  ic) .  This  has  been  optimized  to  minimize  the 
%  time  to  generate  the  3D  ellipse.  The  changes  resulted  in  a 
%  94%  decrease  in  the  cpu  time  required  to  generate  the  ellipse. 

k  =  ic-1; 

iO  =  1;  il  =  diameter  +  1; 
jO  =  1;  jl  =  diameter  +  1; 
while  k  <  diameter+1 
k  =  k  +  1; 

%  This  loops  over  all  the  pixels  in  the  first  plane  to  find 
%  the  pixels  belonging  to  the  ellipse 

for  i  =  iO : il 

for  j  =  j  0  :  j 1 

al  =  [i-ic,  j-ic,  k-ic] ; 
a2  =  [ aA2 ,  bA2,  cA2] ; 
cl  =  al*M; 
c2  =  cl . A2 . /a2 ; 
dist(i,j,k)  =  sum(c2); 

end 

end 

%  If  there  are  no  pixels  belonging  to  the  ellipse  on  this  plane, 
%  then  go  ahead  and  exit  out  of  the  loop  by  setting  k  to  the 
%  final  plane 

if  sum ( sum (dist ( : , : , k) <=1 ) ) ==0 
k  =  diameter  +  1; 

end 

%  If  this  is  not  the  first  or  last  plane,  then  the  program 
%  is  smarter  about  which  pixels  it  checks  for  ellipse  pixels  in 

if  k  ~=  ic  &&  k  ~=  diameter  +  1 
d  =  dist ( : , : , k-1 )  <=  1; 
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e  =  dist  (  :  ,  : , k)  <=  1  ; 

dl  =  dif f ( sum (d, 1 ) ~=0 ) ;  dmin  =  find(dl  ==  1); 
dmax  =  find(dl  ==  -1); 

if  size (dmin, 2)  >  1;  dmin  =  dmin(l);  dmax  =  dmax (2);  end 

el  =  dif f ( sum (e, 1 ) ~=0 ) ;  emin  =  find(el  ==  1); 
emax  =  find(el  ==  -1)  ; 

if  size (emin, 2)  >  1;  emin  =  emin(l);  emax  =  emax (2);  end 

if  dmin  -  emin  <  0; 
j  0  =  emin  -  1 ; 
jl  =  jO  +  3  +  (dmax-dmin); 

if  jl  >  diameter  +1;  jl  =  diameter  +  1;  end 

end 

if  dmax  -  emax  >  0; 
j 1  =  emax  +  1 ; 
jO  =  jl  -  3  -  (dmax-dmin); 
if  jO  <  1;  jO  =  1;  end 

end 

dl  =  dif f (sum (d, 2) ~=0) ;  dmin  =  find(dl  ==  1); 
dmax  =  find(dl  ==  -1); 

if  size (dmin, 1)  >  1;  dmin  =  dmin(l);  dmax  =  dmax (2);  end 

el  =  dif f (sum (e,  2) ~=0)  ;  emin  =  find(el  ==  1); 
emax  =  find(el  ==  -1); 

if  size (emin, 1)  >  1;  emin  =  emin(l);  emax  =  emax (2);  end 

if  dmin  -  emin  <  0; 

10  =  emin  -  1; 
il=i0+3+  (dmax-dmin) ; 

if  il  >  diameter  +1;  il  =  diameter  +  1;  end 

end 

if  dmax  -  emax  >  0; 

11  =  emax  +  1; 
i0=il-3-  (dmax-dmin) ; 
if  iO  <  1;  iO  =  1;  end 

end 

end 

end 


%%  Generate  whole  ellipse 

%  This  section  uses  symmetry  of  the  ellipse  about  the  center  plane 
%  to  generate  the  entire  ellipse  in  matrix  y.  The  defined  distance 
%  matrix  in  the  previous  section  has  a  threshold  of  1.  Distances 
%  less  than  or  equal  to  1  belong  to  the  ellipse;  all  other  distances 
%  do  not . 
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y  =  dist  <=  1; 

x  =  y ( 1 : diameter+1 , 1 : diameter+1 , ic  +  1 : diameter+1 ) ; 
xl ( 1 : diameter+1 , 1 : diameter+1 , 1 : ic-1 )  =  0; 
for  i  =  1 : size (x, 3) 

xl ( : , : , size (x, 3) -i+1)  =  flipud (fliplr (x ( : , : , i) ) ) ; 

end 

y ( 1 : diameter+1 , 1 : diameter+1 , 1 : ic-1 )  =  xl; 

%%  Show  images 

%  For  debug  purposes,  change  the  'while  O'  to  'while  1'  for 
%  of  images  through  the  ellipse. 

while  0 

for  i  =  1: length (y) 

image_view ( 1-y ( : , : , i) ) ; 
pause  (0.15) 

end 

end 
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Intentionally  left  blank. 
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Appendix  C.  Other  MATLAB  Functions 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  image_view ( I ) 


%  Mark  Tschopp 
%  2016 

%  Quick  image  view  function  that  doesn' 
%  the  MATLAB  Image  Processing  toolbox 

imagesc ( I )  ; 
colormap (gray)  ; 
axis  off; 
axis  equal; 
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%  created  by:  Quan  Quach 
%  date:  11/6/07 

%  this  function  erases  any  empty  sheets  in  an  excel  document 
function  deleteEmptyExcelSheets (fileName) 

%  the  input  fileName  is  the  entire  path  of  the  file 

%  e.g.,  fileName  =  'C:\Documents  and  Settings\matlab\myExcelFile.xls' 

excelObj  =  actxserver (' Excel . Application ')  ; 

%  opens  up  an  excel  object 

excelWorkbook  =  excelObj .workbooks .Open (fileName) ; 
worksheets  =  excelObj . sheets ; 

%  total  number  of  sheets  in  workbook 
numSheets  =  worksheets . Count ; 

count=l ; 

for  x=l : numSheets 

%  stores  the  current  number  of  sheets  in  the  workbook 
%  this  number  will  change  if  sheets  are  deleted 
temp  =  worksheets . count ; 

%  if  there's  only  one  sheet  left,  we  must  leave  it  or  else 
%  there  will  be  an  error, 
if  (temp  ==  1) 
break; 

end 

%  this  command  will  only  delete  the  sheet  if  it  is  empty 
worksheets . Item (count) .Delete; 

%  if  a  sheet  was  not  deleted,  we  move  on  to  the  next  one 
%  by  incrementing  the  count  variable 
if  (temp  ==  worksheets . count) 
count  =  count  +  1; 

end 

end 

excelWorkbook . Save; 
excelWorkbook .Close (false)  ; 
excelObj . Quit ; 
delete (excelObj ) ; 
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